#!/bin/csh
# 2d_timeprofiles.csh
# Plot the time evolution of the horizontally averaged profile of a given field 
# The script needs the file time_profiles.dat in the output directory, which 
# must be produced using make_timeprofiles.csh 
#
# Run it using ./2d_timeprofiles.csh field
# where field can be either temperature, viscosity, composition, 
# expansivity or conductivity
###################################################################

source simulation.dir

if (! -d $dirout) then
  echo $dirout 'does not exist. Check the path stored in the variable "dirout" '
  exit
endif

set dirplots = {$dirout}Plots/ 

set makejpg = 1
set jpgdensity = 150
set deletepdf = 1

#------------------------------------------------------
# Read the field to plot from standard input
#------------------------------------------------------
set field = ${1}

set field = ${1}
if (${1} !~ temperature && ${1} !~ composition && ${1} !~ viscosity && ${1} !~ conductivity && ${1} !~ expansivity) then
  echo 'The field you want to plot is either not (yet) supported in this routine or has been misspelled'
  exit
endif

#------------------------------------------------------
# Create plots directories if necessary
#------------------------------------------------------
if (! -d $dirplots) then
  mkdir $dirplots
endif
if (! -d $dirplots$field) then
  mkdir $dirplots$field
endif

#####################################
# GMT plot
#####################################

#----------------------------------
# Set gmtdefaults
# ---------------------------------
source setgmtdefaults

# Select the input file obtained from make_timeprofiles.csh
set tzfile = {$dirout}time_profiles.dat
echo "Output file:"
echo $tzfile

set X0 = 4
set Y0 = 5

set dimx = 16
set dimz = 8
set cscalex = `echo "$dimx/2" | bc -l`
set cscalez = -2
set cscaled = `echo "$dimx*0.8" | bc -l`

set xlabel = "Time"
set zlabel = "Depth"

# Depth range
set zmin = 0
set zmax = 1

# Determine max and min time
set time0 = `awk 'NR > 1 {if(min==""){min=$1}; if($1< min) {min=$1}} END {printf "%12.5f", min}' $tzfile`
set time1 = `awk 'NR > 1 {if(min==""){min=max=$1}; if($1>max) {max=$1}} END {printf "%12.5f", max}' $tzfile`
set atime = `echo "scale=3; ($time1 - $time0)/5" | bc -l`

# time and z increments for grid contours
set dt = `echo "scale=5; ($time1 - $time0)/500" | bc -l`
set dz = 0.002

#----------------
# Temperature 
#----------------
if ($field == temperature) then
  echo "Plot temperature" 

  set output = temperature.ps
  set pdfoutput = temperature_tztime.pdf
  set jpgoutput = temperature_tztime.jpg

  set clabel = "Temperature"
  set cscale = polar
  set ctick = 0.5
  set stepcontour = 0.02
  set opt = ""

  set min = 0
  set max = 1
  set delta = `echo "scale=2; ($max - $min)/10" | bc -l`

  awk 'NR>1 {print $1, 1-$2, $4}' $tzfile >! p

endif

#----------------
# Composition
#----------------
if ($field == composition) then
  echo "Plot composition"

  set output = composition.ps
  set pdfoutput = composition_tztime.pdf
  set jpgoutput = composition_tztime.jpg

  set clabel = "Composition"
  set cscale = wysiwyg
  set stepcontour = 0.05
  set ctick = 0.5
  set opt = "-I"
  set adaptive_cscale = 0

  set min = 0
  set max = 1
  set delta = `echo "scale=2; ($max - $min)/10" | bc -l`

  awk 'NR>1 {print $1, 1-$2, $7}' $tzfile >! p

endif

#----------------
# Viscosity
#----------------
if ($field == viscosity) then
  echo "Plot viscosity" 

  set output = viscosity.ps
  set pdfoutput = viscosity_tztime.pdf
  set jpgoutput = viscosity_tztime.jpg

  set clabel = `echo log@-10@-Viscosity`
  set cscale = wysiwyg
  set opt = "-I"

  set min = `awk 'NR > 1 {if(min==""){min=log($5)/log(10)}; if(log($5)/log(10)< min) {min=log($5)/log(10)}} END {printf "%12.5f", min}' $tzfile`
  set max = `awk 'NR > 1 {if(min==""){min=max=log($5)/log(10)}; if(log($5)/log(10)>max) {max=log($5)/log(10)}} END {printf "%12.5f", max}' $tzfile`
  set delta1 = `echo "scale=1; ($max - $min)/20" | bc -l`
  set delta2 = `echo "scale=1; ($max - $min)/5" | bc -l`
  set delta = `echo "scale=1; ($max - $min)/20" | bc -l`

  set stepcontour = $delta1
  set ctick = $delta2

  awk 'NR>1 {print $1, 1-$2, log($5)/log(10)}' $tzfile >! p

endif

#----------------
# Expansivity
#----------------
if ($field == expansivity) then
  echo "Plot expansivity" 

  set output = expansivity.ps
  set pdfoutput = expansivity_tztime.pdf
  set jpgoutput = expansivity_tztime.jpg

  set clabel = `echo Thermal expansivity`
  set cscale = wysiwyg
  set opt = "-I"

  set min = `awk 'NR > 1 {if(min==""){min=$11}; if($11< min) {min=$11}} END {printf "%12.5f", min}' $tzfile`
  set max = `awk 'NR > 1 {if(min==""){min=max=$11}; if($11>max) {max=$11}} END {printf "%12.5f", max}' $tzfile`
  set delta1 = `echo "scale=2; ($max - $min)/20" | bc -l`
  set delta2 = `echo "scale=2; ($max - $min)/5" | bc -l`
  set delta = `echo "scale=2; ($max - $min)/20" | bc -l`

  set stepcontour = $delta1
  set ctick = $delta2

  awk 'NR>1 {print $1, 1-$2, $11}' $tzfile >! p

endif

#----------------
# Conductivity
#----------------
if ($field == conductivity) then
  echo "Plot conductivity" 

  set output = conductivity.ps
  set pdfoutput = conductivity_tztime.pdf
  set jpgoutput = conductivity_tztime.jpg

  set clabel = `echo Thermal conductivity`
  set cscale = wysiwyg
  set opt = "-I"

  set min = `awk 'NR > 1 {if(min==""){min=$12}; if($12< min) {min=$12}} END {printf "%12.5f", min}' $tzfile`
  set max = `awk 'NR > 1 {if(min==""){min=max=$12}; if($12>max) {max=$12}} END {printf "%12.5f", max}' $tzfile`
  set delta1 = `echo "scale=2; ($max - $min)/20" | bc -l`
  set delta2 = `echo "scale=2; ($max - $min)/5" | bc -l`
  set delta = `echo "scale=2; ($max - $min)/20" | bc -l`

  set stepcontour = $delta1
  set ctick = $delta2

  awk 'NR>1 {print $1, 1-$2, $12}' $tzfile >! p

endif

surface p -Gp.grd -I$dt/$dz -R${time0}/${time1}/$zmin/$zmax
makecpt -C$cscale -T$min/$max/$stepcontour $opt -Z > p.cpt
 
grdimage p.grd -X$X0 -Y$Y0 -R -JX${dimx}/-${dimz} -Cp.cpt -Ba{$atime}:${xlabel}:/a0.2:${zlabel}:WSne:."": -K >! $output
grdcontour p.grd -JX -C${delta} -W0.6/50/50/50 -O -K >> $output

psscale -D$cscalex/$cscalez/$cscaled/0.3h -Cp.cpt -N -Ba${ctick}:"${clabel}": -O >> $output

#----------------------------------------
# Conversions and cleaning
#----------------------------------------

ps2raster $output -Au -P -Tf -F$pdfoutput
if ($makejpg == 1) then
  ps2raster $output -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  mv -f *.jpg $dirplots$field
endif

if ($makejpg == 1 && $deletepdf == 1) then
  rm -f *.pdf
else 
  mv -f *.pdf $dirplots$field
endif

rm -f *.ps *.grd *.cpt p* 

